Collective oscillations in driven coagulation 
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We present a novel form of collective oscillatory behavior in the kinetics of irreversible coagulation 
with a constant input of monomers and removal of large clusters. For a broad class of collision rates, 
this system reaches a non-equilibrium stationary state at large times and the cluster size distribution 
tends to a universal form characterised by a constant flux of mass through the space of cluster sizes. 
Universality, in this context, means that the stationary state becomes independent of the cut-off 
as the cut-off grows. This universality is lost, however, if the aggregation rate between large and 
small clusters increases sufficiently steeply as a function of cluster sizes. We identify a transition to 
a regime in which the stationary state vanishes as the cut-off grows. This non-universal stationary 
state becomes unstable, however, as the cut-off is increased and undergoes a Hopf bifurcation. After 
this bifurcation, the stationary kinetics are replaced by persistent and periodic collective oscillations. 
These oscillations carry pulses of mass through the space of cluster sizes. As a result, the average 
mass flux remains constant. Furthermore, universality is partially restored in the sense that the 
scaling of the period and amplitude of oscillation is inherited from the dynamical scaling exponents 
of the universal regime. The implications of this new type of long-time asymptotic behaviour for 
other driven non-equilibrium systems are discussed. 
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PACS numbers: 82.40.Bj,82.40.Ck,83.80. Jx 

The statistical dynamics of irreversible coagulation 
have been studied for almost a century since the pio- 
neering work of Smoluchowski on Brownian coagulation 
of spherical droplets. See [1] for a modern review. It 
nevertheless remains an important branch of statistical 
physics. This is in part due to its status as a paradigm 
of non-equilibrium kinetics, but primarily due to its con- 
nections to variety of important modern problems. We 
particularly highlight applications in cloud physics 
surface growth 01 and planetary physics y^. In these 
examples, coagulation of clusters is supplemented with a 
source (or effective source in the case of of small clus- 
ters or "monomers" . Such driven coagulation, in which 
monomers are supplied to the system at a constant rate, 
is the main focus of this article. One may expect the 
kinetics of such a system to become stationary for large 
times jSi] with the loss of clusters due to coagulation com- 
pensated by the supply of new clusters provided by the 
input of monomers. We show below that this intuitive 
picture is not always correct and demonstrate the possi- 
bility of a new and strikingly different long time behavior 
characterised by time-periodic oscillatory kinetics. 

Before we begin, let us introduce a large mass cut-off, 
M. Above this size clusters are removed from the system. 
Physically this could be literal removal as in the case of 
large droplets preferentially precipitating out of a cloud. 
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or quenching of reactivity due, for example, to charge 
accumulation. Our primary motivation for introducing 
it, however, is theoretical and we shall focus on what 
happens as M — > oo. The basic quantity of interest is 
the cluster size distribution denoted by Nm{t). It is the 
average density of clusters of mass m at time t. Assuming 
that the system is statistically homogeneous, N.m{t) has 
no spatial dependence. We denote the coagulation rate 
between clusters (or coagulation "kernel") by X (mi, 7712). 
Suppressing the i-dependence of Nm{t) for brevity, the 
mean-field kinetics satisfy Smoluchowski's equation: 



dtNr, 
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dmiK{mi, m — mi)NmiN„ 
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dmiK{m, mi)Nmi + J S{m — 1) 



Dm [Nn 



where 



/■M 

Dm [N,n] =N^ dmi K{m, mi) N^, (2) 

J M-m 

removes clusters larger than M and J is the monomer 
injection rate. We study the family of kernels 

K(mi,m2) = ^ (mim^ -|- m'^m^) , (3) 

which includes many of the commonly studied models [ij . 
Eq. ([3]) can also capture the asymptotics of most physi- 
cally relevant kernels. We mostly consider cases for which 
^ + v < 1 . This avoids complications due to gelation [l[ . 
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The stationary solution of Eq. ([T]) without cut-off was 
found in It is a power law for large m: 



N„ 



J[l-{v- ^)2] cos[7r(i/ - ^)/2] 
47r 



TO 2 . (4) 



The exponent 



implies a constant flux of mass 



through the space of sizes, to. It is a standard example 
of a non-equilibrium stationary state with a conserved 
current. From Eq. (U), this stationary state exists only 
if |:/ — ul < 1, a fact which is true for any scale invariant 
kernel |7|. One might ask what happens if — /i| > 1? 
This can occur in practice. Examples include coagula- 
tion of ice clusters in planetary rings ^^l, gravitational 
clustering jSj] and droplet sedimentation in static fluids 



The fact that the constant flux stationary state only 
exists for a certain class of kernels has long been appre- 
ciated in the theory of wave kinetics There, the 
constraint \v — ^\ < 1 would be interpreted in terms 
of universality. If one solves the stationary version of 
Eq. ([1]) with finite cut-off, M, and studies the behavior 
as M — > oo one finds that when — < 1, the leading 
order terms becomes independent of M as M ^> oo. The 
stationary state thus tends to the above universal form 
found in [6j. If, on the other hand, \v — fi\ > 1, the sta- 
tionary state is non-universal and retains a dependence 
on M as M — > oo. This phenomenon is referred to as 
nonlocality of interaction (in the mass space) in the sense 
that all masses remain strongly coupled to the largest 
and smallest masses in the system. By extension, the 
interactions in the regime \v — ii\ < 1 are termed local 
although this is a rather weak form of locality. The pres- 
ence of a finite cut-off is essential to obtain a stationary 
state in the nonlocal regime as discussed in 



Almost nothing is presently known about the shape of 
Njn in the nonlocal regime. We developed an algorithm 
to compute the exact stationary solution of the discrete 
version of Eq. ^ with cut-off by converting it into a two- 
dimensional minimisation problem which can be easily 
solved numerically for modest values of M. For details 
see Appendix |^ Some typical results are shown by the 
symbols in Fig. [TJ It is clear that the nonlocal stationary 
state is not a simple power law. To obtain some analytic 
understanding, one possible way forward was outlined in 
[I2]. If clusters of size to grow primarily by interaction 
with clusters of mass mi <C to, which is the essential 
feature of nonlocal interactions, one can Taylor expand 
the righthand side of Eq . ([Ij and obtain an almost linear 
equation for Nmit) [12j]. The dominant terms in this 
equation are 



dN„. 

dt 



(5) 




FIG. 1. Comparison of asymptotic approximation, Eq. ([6]) 
(solid lines) to the true stationary state of Eq. IT} (symbols) 
with the kernel given by Eq. ^ for several values of 1/ and fi 
chosen in the nonlocal regime. The cut-off is M = 10*. 



where the t dependence of Nm has been suppressed and 



= 1^ 



Nrmdmi- 



M 



miNjn-^dmi. 



Extension of the limits of integration of these latter in- 
tegrals to M and 1 respectively is a further assump- 
tion which needs to be justified a-posteriori. The self- 
consistent calculation detailed in '13] for the case fj, — is 
easily extended to obtain the following stationary asymp- 
totic solution of Eq. ^ in the limit of large M: 



TV,; - v/27 J log {M)M~^Ivr ''m' 



(6) 



where ^ ~ v — adopting the convention that v > 

in Eq. (jSj. Detailed derivations of Eqs. ([5]) and dH) are 
provided in Appendices [B] and [C] Equation ([5]) approxi- 
mates well the true stationary state as indicated by the 
solid lines in Fig. [T] Note that there are no adjustable 
parameters. A striking feature of Eq. ([6]) is that the 
prefactor of the stationary state vanishes as M 00 re- 
flecting the non-universality of the nonlocal regime. Sim- 
ilar behaviour was observed in the instantaneous gelation 
regime in [isj although there is no gelation here. 

The vanishing of the stationary state in the limit 
M ^ 00 poses a conceptual problem since it suggests 
the removal of the conduit linking the source to the sink. 
In order to investigate how the conserved mass current 
is carried in the nonlocal regime, we computed dynami- 
cal solutions of Eq. ([T]) in the nonlocal regime using the 
numerical algorithm developed in [T3]. The results were 
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FIG. 2. Main panel: Total mass vs time for different values 
of M with V = —ix= |. Inset: Collapse obtained by rescaling 
the data according to Eqs. (O. 



surprising. For small values of M, the numerical solu- 
tion converged to the exact stationary state as expected. 
Once M exceeded a certain value, however, the numerical 
solution never reached the stationary state. The typical 
behaviour of the total mass as a function of time for dif- 
ferent values of M is shown in the main panel of Fig. [2] 
for the case v = —fi = |. Stationarity is reached only 
for smaller values of M. For larger M we observe col- 
lective oscillations which seem to persist indefinitely (we 
stopped the computation after several hundred periods) . 
The period and amplitude grow with M . 

The intriguing possibility thus arises that the station- 
ary state becomes unstable as M increases. Our algo- 
rithm for computing the stationary state is not dynami- 
cal and makes no distinction between stable and unsta- 
ble fixed points. We therefore input the exact stationary 
state as an initial condition for the dynamical code and 
added a small perturbation. The results for the density 
are shown in the inset of Fig. [3] The perturbation grows 
to a finite amplitude in a clear indication of instability. 
A lin-log plot of the amplitude of the successive max- 
ima of the perturbation, as shown in the main panel of 
Fig.[3l indicates exponential growth, a clear sign of linear 
instability. We used Mathematica to compute the eigen- 
value, Cmax, of the linearization of the discrete version 
of Eq. ([T]) about the stationary state having maximum 
real part. This analysis confirmed the instability. The 
growth rate agrees well with numerics (see main panel of 
Fig. [3]). For fixed ly and /i, the stationary state under- 
goes a Hopf bifurcation as M is increased. The eigenvalue 
Cmax crosses the imaginary axis at a critical value of M 
(see inset of Fig. S]) giving birth to a limit cycle and os- 
cillatory behavior. The structure of the instability as a 
function of v and fi for fixed M is non-trivial as shown 
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FIG. 3. Numerical evolution of a perturbation of the sta- 
tionary state for v = —fi = | and M = 100. Main panel: 
amplitude of successive maxima of the perturbation (circles). 
The solid line is the prediction of linear stability analysis. 
Inset: oscillations of the total density. 
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FIG. 4. Main panel: Re[^]inax for kernels fx = —v plotted as 
a function of u for different values of the cut-off, M. Inset: 
R.e[(^]niax as a function of M for different values of i/. 

in the main panel of Fig. SI For fixed M, the stationary 
state becomes stable again for sufficiently large values 
of i^, a fact for which we have no intuitive explanation 
at present. Such limit cycles appearing in mean-field 
equations can be destroyed by noise [15]. To check the 
robustness of this phenenomenon, we performed Monte- 
Carlo simulations of the Markus-Lushnikov model (see 
p^ ) with a source and sink of particles. Typical results 
are shown in Fig. [S] Oscillations are clearly visible which 
remain coherent in the presence of noise. 
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FIG. 5. Total mass vs time in a Monte-Carlo simulation of the 
Markus-Lushnikov model with a source and sink and kernel 
given by Eq. © with ^= -v ^ -0.95 and M = 300. 



To understand nonlinear aspects of the instability such 
as the period and amplitude of nonlinear oscillation we 
return to Eq. ([T]). Each period corresponds to a pulse of 
mass through the space of sizes. A movie provided with 
the arxiv version of this paper illustrates these pulses. 
For details of the parameters see Appendix |D] Each pulse 
almost resets the mass of the system to zero as evident 
from the main panel of Fig. [21 Let us suppose each pulse 
grows with self-similar size distribution, 



(7) 



where s{t) is a typical size and a is an exponent to be 
determined. Substituting Eq. ([7]) into Eq. ^ and bal- 
ancing dependences on t requires that s = s'^+/'+o+2 
Since the mass contained in each pulse grows linearly in 



time, mNm{t)dm = 
differentiating gives s s 



Jt. Substituting Eq. (O and 
Consistency requires 



The period is estimated as the time, tm, required for the 
typical mass, s{t), to reach M. The amplitude. Am, is 
estimated as the mass supplied in one period. We thus 
obtain the following scalings for tm and A]\j with M: 



TM ~ M- 



A 



M 



J M- 



(9) 



These scalings are verified by the data collapse presented 
in the inset of Fig. [5J Universality is in a sense restored 
since the earlier universal behavior of Eq.Q can now 
be understood as the special case in which F{£) has the 
special form which cancels s{t) from Nm{t) in Eq.Q. 

We believe that the phenomena presented here are un- 
likely to be restricted to coagulation. Many driven dis- 
sipative systems with conserved currents must satisfy a 
locality criterion analogous to the one discussed here 17 1 
and may be candidates for oscillatory behaviour when 
this criterion is violated. In particular, the kinetic equa- 
tion for isotropic 3- wave turbulence, which is closely anal- 
ogous to Eq. (QJ, becomes nonlocal when — /i| > 3 iisj. 
Furthermore, the oscillatory behaviour discussed in this 
article may even have been already observed experimen- 
tally in measurements of non-equilibrium phase separa- 
tion of binary mixtures with slowly ramped temperature 
[20j . In this system, droplets of one phase coagulate 
inside another during demixing with nucleation provid- 
ing the source of "monomers" although the coagulation 
process is not obviously nonlocal in our sense. This nev- 
ertheless seems like a potentially fruitful direction for fur- 
ther investigation since the theory presented here makes 
several testable predictions about the oscillatory kinetics. 



sit) ~ t^^^. 



(8) 
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Appendix A: Algorithm for finding exact stationary 
solution of the Smoluchowski equation 

Consider the discrete form of the stationary Smolu- 
chowski equation: 



^ m — 1 

= - ^ K{mi,m - mi) N„i^ Nm-mi 

mi — 1 
M 



Nm ^ K{m,mi) 



(Al) 



We use the kernel considered in the main text: 



K{mi, 777,2) — 2 ('T*i"^2 + "^i"^2 ) ■ (-'^2) 

Here g is a constant which can be helpful to keep track of 
dimensions but which is usually set equal to one. Denote 
the p"^ moment of the size distribution by 



M 



mi — 1 



If the moments Aif^ and Aii, were known, we could find 
the solution of (|A1[) by iteration: 



Given that the solution can be expressed in terms 
of the two moments and M,j, the task is now to 
self-consistently determine the values of these moments. 
We can approach this task as a simple two-dimensional 
optimization problem. Eq. (|A3p expresses Nm as a func- 
tion of the pair of moments iVm(A^^, Ai,^). We create an 
objective function, 5'(A^^, A4^), as follows: 



M 



m— 1 
M 



(A6) 



The correct values of and A4^, which we de- 

note by A4i_i* and A4^*, can be found by minimising 



{Mf,*,Mu*) = a.Tg min ^{Mf„M^). (A7) 



This can be done with any numerical minimization 
algorithm. We used the Nelder-Mead downhill simplex 
method. The solution thus obtained is exact to within 
computational error since no approximations have been 
made in formulating this procedure. 

We remark that the problem does not have to be for- 
mulated as an optimization problem. One could treat it 
as two-dimensional root-finding problem. Furthermore, 
by summing Ea. (|Aip . one can derive an independent re- 
lationship between Aifj, and A4i, (in the limit M 00) 
which further reduces the problem to a one-dimensional 
root finding problem. We did experiment with some of 
these alternatives but settled on the procedure described 
above as the most numerically stable and reliable ap- 
proach. 



N — a 



where Gm depends only on the densities of clusters with 
masses less than 777, and is given by 



m— 1 

Gin = ^ X! -^("^1; - "il) ^mi N„ 



The starting value is obtained by setting m 1 which 
gives the monomer density 



7Vi = 



2J/g 



Appendix B: Derivation of the nonlocal 
(A3) Smoluchowski equation 

We assume without loss of generality that v > jj.. It is 
the combination 7 = v~ which determines the local- 
ity of the stationary solution of Smoluchowski's equation. 
The nonlocal case corresponds to 7 > 0. We can use the 

(A4) 

differential approximation outlined in [12] to describe the 
stationary state in this regime. The first step is to rewrite 
the Smoluchowski equation in a particular form. Terms 
describing interactions between a reference mass, m, and 
masses less than ^ are gathered together in one group. 
Those describing interactions with masses larger than ^ 
(A5) are gathered together in a second group. Splitting the 
integrals appropriately and performing some manipula- 
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tions we obtain: 

dmi [K{rni,m-Tni)Njn-rni- K{mi,rn)Nrn\N^^ 

Jo 

Nm. I dmiK{m, mi)N,ni + J 5{m — mo). (Bl) 

Consider the first term which accounts for all interac- 
tions between clusters of mass m and those having mass, 
mi < Y- If the cascade is nonlocal, these interactions are 
primarily with those clusters having mi <^ -y, in which 
case the integrand is strongly concentrated in the region 
mi <§; We can then Taylor expand with respect to 
mi and neglect all terms of O(m^) or higher to obtain: 



d_ 

dm 



dmiK{m, mi) N^-^ N„ 



M 



(B2) 



With the kernel given by Eq. (IA2[) . we get the following 
equation for the stationary state: 



- Nm I dmiK{m,mi)Nm^_ 
+ J 5{m — mp). 



: 



dm 



[M<+lm^'N„ 



dm 



[X<+im"^7V„ 



2 J 



- X>m^iV„, - Mlm'^Nm + — 5{m - mo),(B3) 

9 

where M.^ and denote the upper and lower partial 
moments: 



M< = 



M> = 



dmim^Nm^ 



mo 
M 



dmim^Nmi ■ 



(B4) 
(B5) 



The dominant terms in this equation when m 3> mo will 
turn out to be 

= [M<+im''iV„,] - M>m''Nm 

H o(m — mo). 

9 

This statement will have to be justified a-posteriori. In 
order to make further progress let us assume that the 
error made by extending the upper and lower limits of 
integration in the partial moments ■M.'^j^i and AAu to M 
and mo respectively is small. This will also have to be 
justified a-posteriori. The resulting equation is: 



= - — [M^.+im'^Nm] - X.m^iV„ 
dm 

H o(m — mo). 

9 

This is Eq.(5) in the main text. 



(B6) 



Appendix C: Asypmtotic solution of nonlocal 
Smoluchowski equation 



Eq. (|B6P is a linear equation and can be readily inte- 
grated to give 



N„ 



Ce- 



where /3 is a ratio of moments 



M 



(CI) 



(C2) 



and C is a constant of integration. The non-trivial aspect 
of the problem is that the moments AAi, and M.^+i must 
be determined self-consistently from this solution. This 
cannot be done analytically but an asymptotic solution 
for large cutoff, M, can be found which we now describe. 
A general moment of order a is 



M 



Mc^^C I dm m"-''e-"' 

/mo 



t-'^-^^e*dt (C3) 



M-1 



where we have introduced the shorthand notation for 
the combination 



1 



+ a 



7 



Since fi > v + 1 m the nonlocal regime, we would expect 
the moment to grow faster than A^^+i as the cutoff, 
M, is increased. We therefore expect /3 to grow as M 
grows. Let us suppose that it does not grow faster than 
. If this is the case then the upper limit of the integral 
in Eq. (jC3|) tends to infinity as AI grows while the lower 
limit tends to zero. We are therefore interested in the 
behaviour of the integral 

.A 

/(e,A,Ca) / r^-'^e'dt 



as e ^ and A — i> cxd. This integral clearly diverges at 
its upper limit regardless of the value of Ca ■ The leading 
order behaviour as the upper limit grows is 



as A 



(C4) 



It is divergent at its lower limit if Cq > 0. The leading 
order behaviour as the lower limit goes to zero is 



as e 



0. 



(C5) 



The moments of immediate interest correspond to a = 
and a = ^ + 1 which give values for Ca of ^ and 
respectively. For a — e always get a divergence at the 
lower limit. For a = -1-1 we get a divergence at the lower 
limit for 7 in the range < 7 < 1. With this knowledge 



7 



in mind, our task is now to substitute Ea. (|C3l) into the 
consistency condition, Eq. (|C2p and attempt to balance 
the divergences as M (and thus (3) tends to infinity. One 
solution is to balance the divergence coming from the 
lower limit of A^^+i with the one coming from the upper 
limit of A^^. This gives, after some work 



p ^ 7 mn in 

mo 



as M — oo. 



(C6) 



which is consistent with our assumption that /3 should 
grow with M but slower than Iv'P . We can now substi- 
tute this value for /? into Eq. (jC3P and obtain the leading 
order behaviour of Mv and Mfi+i. We find that Mu is 
dominated by its upper limit and grows as 



My Cnio ( — 
mo 



as Ad — > oo. 



(C7) 



On the other hand A^^+i is dominated by its lower limit 
and grows as 



M 



mo 



mo 



as M oo. 
(C8) 

These estimates justify our replacement of the partial 
moments with full moments in the derivation of Eq. (jB6p . 

It remains to find the constant C. This can be done 
by requiring that the total mass flux leaving the system 
is equal to the input flux, J: 

I'M 

J— dmra I dmiK{m,mi) NmNmi- (C9) 

J niQ JM — ra 

Substituting the kernel Ea. (IA2p into this gives two terms: 
= / dm m'^'^'^ N,n / dmim'( N^^ 

9 Jmo J M-m 

I'M I'M 

-t- / dmm'^^^ N,n / dmim'^ Nmi- 

J mo M —m 



With some further analysis one finds that the second 
term is much smaller than the first term as M grows. 
Extending the regions of integration of the partial mo- 
ments as before we obtain the estimate; 



2 J 

— ^ Mn+i My as M oo. 
9 



Using Eqs. ([C7)) and ((C8|) we obtain 



mo 



C = — — m2 ^ In I — 



(CIO) 



Putting this together with Eqs. (jCip and (IC6p we finally 
obtain: 



N* 




(Cll) 



\mo/ V"^o, 
The explicit dependence on the monomer mass, mo 

(which we usually take equal to 1) has been retained in 

order to make the dimensional correctness of the formula 

clear. Setting mp = 1 gives Eq. (6) in the main text. 



Appendix D: Comment on the accompanying movie 

The movie accompanying the arxiv version of this pa- 
per shows the time evolution of the density contrast, 
N{m, t) /N* (m) , relative to the stationary state as a func- 
tion of cluster size, m, in the oscillatory regime. This 
quantity would be 1 if the stationary state were sta- 
ble. Both axes are linear. The movie was generated 
by solving the discrete Smoluchowski equation (without 
coarsegraining) with ly — —fj, = |, a monomer input rate 
of J = 1 and a cut-off of M = 100. 



